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A simple procedure is presented which leads to a dramatic simplification in 
Q_i! the calculation of the relic density of stable particles in the Universe. 
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Any stable species \ contributes to the total mass-energy density in the Universe. 
If its number density cannot be reduced efficiently enough before it decouples from 
the thermal equilibrium, its relic abundance ^ x h^ can be sizeable and it can affect 
the evolution and the age of the Universe. A conservative estimate that the Universe 
is at least 10 billion years old requires Q x hl < JItot^o < 1 0- Furthermore, there is 
growing evidence for dark matter at both galactic and larger length scales Q which 
would most likely require the existence of some type of exotic neutral particles. Since 
such particles are often present in many theories beyond the Standard Model, it is 
important to develop a simple practical procedure for calculating their relic abundance 
with enough precision. 

Two groups |fj, ||] have developed equivalent frameworks for properly calculating 
the relic density of x's, including relativistic corrections. The method of Ref. [[J is 
in practice applicable away from poles and new final-state thresholds, which is most 
often the case. In Ref. also the vicinity of poles and thresholds has been carefully 
studied. Essentially, one needs to calculate the thermally averaged product of the \X 
annihilation cross section and their relative velocity (crv Te \) . In practice, one expands 
(av Te i) = a + bx + 0(x 2 ) in powers of x = T/m x = 0(1/20) in order to avoid difficult 
numerical integrations and approximates (<7v Te \) by a and b. Both techniques give 
equivalent results H |j] in the overlapping region (away from poles and thresholds). 

Unfortunately, in practice the actual calculation of even the first two terms of 
the expansion is typically very complicated and tedious. In this Letter I report on a 
dramatic technical simplification in practical applications of the method of Ref. ||. 

Consider an annihilation of particles x, X m t° a two-body final state. Furthermore, 
in many cases of interest, the final state particles have equal mass mp. (A general 
case of unequal final state masses will be presented elsewhere ||.) Let the momenta 
of the two initial states Xi X-, an d the two final states F, F, be pi, p 2 and ki, k 2 , 
respectively. Srednicki, et al, introduce the function w(s) defined as 



wis 



- J dLISP\M\ 2 = E 1 E 2 av xeh (1) 



where d LISP in this case takes the form 

.lisp = ^ + n-k,- m 4?' Ei( J 2 V (2) 

where Ei )2 = k\ 2 + and | M. | 2 is the square of the reduced matrix element for the 
annihilation process XX ~^ FF summed over the spins of the final-state particles and 
averaged over the spins of the initial particles 0. 
Denote 

/ = \M\ 2 - (3) 
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In general / = f(pi,P2, k\, k?). The integral Eq. ([[]) can be conveniently evaluated 
in the centre-of-mass frame (CM) in which f> 2 = —p\ and k 2 = —k\. After a few 
elementary steps one obtains 



w(s) 



2 6 



1- 
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dcos9f(k2 = —ki, \ki\ = Js/4 — m 2 F ). 



Using the kinematic relation between the Mandelstam variables 

t = (pi - h) 2 

one can express Eq. (|) as 
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w(s) 
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t+(s) 



2 5 vr Sv /l -4m2/iA-W 
where £± = t(cos6 = ±1) = to =t ^> 



dtf(s,t), 



t a (s) 



m x + m F -- 



(6) 



(7) 



and 



At(s) 



s 



'1- 



4m^ 



s V s 

Notice that one can always express / as a function of the Mandelstam variables s and t 
only ||. In particular, u can be eliminated by using the relation s+t+u = 2m^+2m F . 
Furthermore, in calculating the relic density it is convenient to introduce |§ 

s 



z 



4m^ 



in terms of which Eq. (Bf) can be simply rewritten as 



w(z) 



1 



t+(z) 



dtf(z,t). 



(9) 



(10) 



2 7 Txm 2 x z\f\ — z Jt-(z) 

In order to calculate the relic abundance of x's one needs to solve the Boltzmann 
(rate) equation. The actual quantity that appears in the rate equation is the thermally 
averaged product (<rv Te {), which is usually approximated by a+bx, as mentioned above. 
One of the main results of Srednicki, et al. 0, was to show that 



— w(z = 1) 



b = — 



3 dw(z 
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dz 



3w(z = 1) 



(11) 
(12) 
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While these formulae look deceptively simple, the actual calculations can be, and 
in practice are, very cumbersome and often virtually unmanageable. For example, 
in a relatively simple case of calculating the interference term between two t-channel 
amplitudes (with the masses of the exchanged particles denoted by /ii and ^2) one 
needs to evaluate several integrals of the type //+ dtt n (t — /if) _1 (t — /uf) -1 ' an d i R 
general in calculating w(z) one has to compute a whole multitude of more complicated 
integrals. Once such a (lengthy) expression for w(z) is found, one needs to next 
take a derivative dw/dz which in general leads to even more cumbersome formulae. 
Additional highly non-trivial computational complications arise when the masses of 
the final state particles are not equal. 

Below I show that a great deal of these difficulties can be avoided. In fact, one 
can avoid performing any integrals completely. First, one can always conveniently 
express || a and b in terms of "reduced" variables a and 6° 



yjl — nip I w? x a 
yjl — nip I m 



-3 + 



2 / 2 
m F /m 



a° + b° 



4 (1 — nip/m x ) 

I will show below that a and 6° can be written in a simple and elegant form as 
1 



(13) 
(14) 
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Tcm: 



-3m 



df(z = 1) 



- v + 2 K _ 2 ^/(* = i) , 3 0/(* = i; 

X. * A 



dt 



m, 



dt 2 



+ 



dz 



(15) 
(16) 



and f(z = 1) should be understood as f(z = l,t (z = 1)), etc. Equations (fL5|) 
and (0) are the main result of this Letter. They can also be readily generalized to 
the case of unequal final-state masses ||. Higher order terms of the expansion can 
also be easily derived. It is clear that the whole procedure of calculating a and 6° 
has now been reduced to merely writing down \M.\ 2 , substituting all the variables 
in |7W| 2 in terms of z and t and next taking a few relatively simple derivatives. In 
fact, one can easily do all these steps entirely with the help of any advanced algebraic 
program. The truly difficult part of computing a and 6° - performing complicated 
integrations in deriving w(z) - has been completely eliminated. 

In order to prove Eqs. (|15"D and (0) notice that for any regular function g(z, t(z)) 
and its integrand G(z) = J dt g(z,t(z)) one can show that 



lim 

Z->1 



1 



lim 

Z->1 



— z Jt-(z) 



t+(z) 



dtg(z,t(z)) 



Am 2 x ^l-m F /mlg(z = 1, t (z = 1)) (17) 
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where I have used lim^i t± = t Q (z = 1), lirn^i At = 2m 2 yl — m^/m 2 a/1 — z and 
UmG(z,t±) = lim[ G(z, t ) ± g(z, t )At + - dg ^ ^ (At) 2 (18) 



Z-s-1 Z->1 



2 dt 



6 at 2 

(In deriving Eq. (|i"6|) one needs to apply the formula fllTD twice (for / and df/dz), and 
keep terms up to (At) 3 only in the expansion of / in order to cancel the term (z— 1)~ 3 / 2 
which appears in dw/dz in Eq. (0).) Applying these formulae to Eqs. (0), Q13D , 
and flT4]) one obtains Eqs. (|I5|) and (pT6|) after a few simple steps. 

As an example, consider the process XX ~^ ZZ in the Minimal Supersymmetric 
Standard Model (MSSM) with x being the lightest neutralino. The two s-channel 
diagrams include the exchange of the scalars h and H, while the four neutralinos Xi 
(i = 1, 4) (including x = Xi) are exchanged in the t- and ^-channels. Since neutralinos 
are Majorana particles, Eq. (^) now reads 

/ = \ E ElMixx^ZZ)^ 1 - £ E\Ms-M t + M u \ 2 (19) 

helicity spin helicity spin 

— f( ss ) _ f(«t) _|_ f{su) _|_ jr(tt) _ j?(tu) _|_ 

After expressing the external momenta in terms of t, u, and z one finds, e.g., 

helicity spin 

= ^ £ (Off) 2 (Og) ^ _ 1 _ ££tf>M, (20) 

COS V-w ij=i [ z m x o)[t m x a) k =0l=0 



fW(z,t) = 2»(i E ElM^ul 



helicity spin 

= -4" E (OffW ) 2 tt r4 2~T E E /iSfW- (2i) 

cos 4 # w ^ (t - m 2 o ) (w - m 2 o ) ^ ^ 

It is also clear that f( uu \z,u) = f( tt \z,t). In the convention used here O"^ = 
-±(N i3 N j3 - NuNjt) (in the basis (B°, W$, H b , H t )), the matrix N {j is real and the 
neutralino masses can be either positive or negative. The expressions for the coef- 
ficients f£P and fki, introduced above, are rather lengthy and will be given else- 
where 0. I will neglect here also the contribution from the s-channel exchange ||. 
Finally, one obtains 

Q _g\ml-m%) 4 (C^) 2 (C>if 2 



a 



E V V \ " (22) 



4ttcos 4 # w ^ A x oA x o 
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and 



22 — .o .o J x 



64tt cos 4 e w m% A 3 A 3 

A v0 A v0 (48m Y om y omt + 32m Y om v + 32m Y om Y + 32m Y 

AjAj\ A j A j A A j A A, j A A 

00 QO QO /10 

-48m v om v om v m 7 — 16m v om v m 7 — 64m v om v m 7 + 32m^m 7 

Aj Aj A A£ A Aj A A ^ 

+36m v om v om^ + llm v om v m* + 29m x om x m| + \2vrcjm 7j + 32m? 



Xj 



+2A v oA Y o (4)71,0 

Xj Xj \ Aj 



m v om 
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24m x om x qm x m| + 10m x om x m| 



— 22m x om x m| + 20m x m z + 20m x om x (im^m| + 2m x om x m% 



Xj Xj 



+46m x om x m| + 5m x m| — 12m x om x m| — 24m x om x m| — 2hm x m% 



Xi Xj 



+2A v oA^,o (28m Y om Y om' 



Aj Aj 



24m v om v om^m 2 7 

Aj A? A ^ 



Xj X ^ 



— 26m v om^m 7 + AAm^m 7 — Am v om v om v m 7 + 82m v om„m 7 

An A K n A.% A.j A. n Aj A ^ 



m x m 2 



+38m x om x m| — I3m x m z — 2Am x om x m z — 12m x o 
+16A v oA v o (4m v om v omJ — 4m Y om v omJm 2 7 — 6m y om[m 

A-j Xj \ Ai A j A Aj A J A " A j A 



31m„m c 



2 
Z 

-J \ ~- ~J " ■ ~- 'J ' '* 

-6m v omIm 2 7 + Am%m 2 7 — 4m v om v »m v ffl7 + 12m v om^mi 

Aj A A " Aj Aj A Aj A ^ 



+12m v omtmi — 5m^mi + Ara v om v ora 2 m% — 6m v o 

AjA^ A ^ AjAjA"^ \ 

—6m x om x m z — 2m x m z + 3m x m 8 z ) 
+32(A x o + A x o) (m x m| - 2m x m| + m 2 x m\) 



Aj A z 



(23) 



m 2 o . 



where Ajj = m z 

These expressions reduce nicely to the Eqs. (3. 19. a) and (3.19.b) of Ref. 
in the limit in which \ becomes an almost pure anti-symmetric higgsino Ha = 
l/y/2 (0,0, -1, 1) (corresponding in the MSSM to M 2 > fx,m z f§). In this limit the 
second lightest neutralino x° is an almost pure symmetric higgsino H s = 1/^2(0,0,1,1) 
and it is almost mass-degenerate with x, |m x x o| ~ /x(l =f e) and m x o ~ — m x . In 
Ref. || the contribution from the two heavier neutralinos, which in this limit are al- 
most pure gauginos, has been neglected. However, as has been pointed out by Drees 
and Nojiri ]7|], they have in this limit non-negligible higgsino components which also 
contribute to the considered process. The expressions (|22[) and (^) are free from this 
problem because they are valid for any type of neutralino and include contributions 
from all four exchanged neutralinos. 

The method presented here has also been tested on a number of other cases (e.g., 
XX ~^ (Z,h, H, A) — > ff), for which simple analytic expressions are available. A 
complete set of relevant expressions for the MSSM case will be presented elsewhere [j5]. 
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Note Added: After this work had been completed and submitted to a journal, I was 
made aware of Ref. || in which the expansion coefficients for av re \ were derived using 
a different method, without specifying a thermal averaging procedure. When these 
expressions are applied to the method of thermal averaging of Ref. jy, which does 
not include relativistic corrections [§, [!|, they become different from the expression 
presented here by a factor b — b x = — |a 0]. 
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